#include "analyzeREADS.h"
#include "legoland.h"
#include "frequency.h"
#include "kmeans.h"
#include <iostream>
#include <fstream>
#include <cctype>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>

/*************************************************************************
* analyzeREADS class
*   This class loads pre-existing libraries and then analyzes reads
*   supplied by the user in either FASTA or FASTQ format. For each
*   read, the probabilities for all of the taxa in the library are
*   calculated for user specified k values. Results are written to file.
**************************************************************************/

using namespace std;
bool file_exist(char * file_name);

/*************************************************************************
* analyzeREADS -  constructor
*   retrieves user's parameters, loads the library and runs the analysis
**************************************************************************/
analyzeREADS::analyzeREADS()
{
    if(!get_parameters())
        cout << "An error occured. Analysis of reads was not conducted. Please refer to previous error messages.";
    else
    {
        load_library();
        analyze();
    }
}

/*************************************************************************
* analyzeREADS -  deconstructor
F*   deletes all dynamically allocated memory class variables
**************************************************************************/
analyzeREADS::~analyzeREADS()
{
    int i;
    delete reference_library_file;
    for(i=0; i<num_gs; i++) delete taxa[i];
    delete taxa;
    delete file_of_files_reads;
    x_axis_file = NULL;
    delete x_axis_file;
    y_axis_file = NULL;
    delete y_axis_file;
    delete output_file_tag;
    for(i=0; i<(k_max-k_min+1); i++)
        delete l[i];
    delete l;
}

/*************************************************************************
* analyze
*   Determines the file format and if the file supplied by the user is
*   a file of reads or a file of files (each of which contains reads).
*   For each k in the user specified range, analysis is performed.
**************************************************************************/
bool analyzeREADS::analyze()
{
    int k_value;
    char line[5000];
    bool fof;
    short int format_file;          //0=FASTA, 1=FASTQ
    format_file=format(fof);
    if(format_file==-1)
    {
        cout << "The reads need to be writted in either FASTA or FASTQ format.\n Analysis not run.\n";
        return false;
    }

    ifstream in;
    if(fof)
    {
        for(k_value=k_min; k_value<=k_max; k_value++)
        {
            in.open(file_of_files_reads);
            if (!in.is_open())
            {
                cout << "could not open file of reads" << endl;
                return false;
            }
            in.getline(line,5000,'\n');
            map_file(format_file,false,line,k_value);               //deletes any prior files by the raw name
            while(in.peek()!=EOF)
            {
                in.getline(line,5000,'\n');
                map_file(format_file,true,line,k_value);
            }
            in.clear();
            in.close();
        }
    }
    else
    {
        for(k_value=k_min; k_value<=k_max; k_value++)
            map_file(format_file,false,file_of_files_reads,k_value);
    }
    return true;
}

/*************************************************************************
* generate_output_files
*   calls functions to write out results - PUBLIC
**************************************************************************/
bool analyzeREADS::generate_output_files()
{
    int num_reads;
    if(!print_summary_statistics(num_reads))
    {
        //this error is printing out
        cout << "Error printing out summary statistics. Check that raw_read_results files exist.\n";
        return false;
    }
    print_probability_statistics(num_reads);
    return true;
}

/*************************************************************************
* print_probability_statistics
*   parameters:
*       -num_reads: number of reads in the file
*   Prints out aggregate probabilities over all k examined
**************************************************************************/
bool analyzeREADS::print_probability_statistics(int & num_reads)
{
    ifstream in;
    int k_value;
    float val;
    char file[500];
    char line[300];
    int i,j;
    float **total_all_k=new float * [num_reads]; for(i=0; i<num_reads; i++) total_all_k[i]=new float [num_gs];
    for(i=0; i<num_reads; i++)
    {
        for(j=0; j<num_gs; j++)
            total_all_k[i][j]=0.0;
    }

    for(k_value=k_min; k_value<=k_max; k_value++)
    {
        sprintf(file,"%s_raw_read_results_k%d.txt",output_file_tag,k_value);
        in.open(file);
        if(in.is_open()) {for(i=0; i<num_reads; i++) delete total_all_k[i]; delete total_all_k; return false;}
        for(i=0; i<num_gs; i++) in.getline(line,300,'\t');
        in.getline(line,300,'\n');

        j=0;
        while(in.peek()!=EOF)
        {
            in.getline(line,300,'\t');
            for(i=0; i<num_gs; i++)
            {
                in>>val;
                total_all_k[j][i]+=val;
            }
            in.getline(line,300,'\n');      //get final CR
            j++;
        }
        in.clear();
        in.close();
    }

    ofstream out;
    sprintf(file,"%s_over_all_k.txt", output_file_tag);
    out.open(file);
    sprintf(file,"%s_raw_read_results_k%d.txt",output_file_tag,k_value-1);
    in.open(file);
    for(j=0; j<num_reads; j++)
    {
        in.getline(line,300,'\t');
        out << line << "\t";
        for(i=0; i<num_gs; i++)
        {
            in>>val;
            out << total_all_k[j][i] << "\t";
        }
        in.getline(line,300,'\n');      //get final CR
        out << endl;
    }
    in.clear();
    in.close();
    out.clear();
    out.close();
    for(i=0; i<num_reads; i++)
        delete total_all_k[i];
    delete total_all_k;
    return true;
}

/*************************************************************************
* print_summary_statistics
*   parameters:
*       -num_reads: number of reads in the file
*   Prints out number of reads mapped and the number of reads mapped
*   to each taxa in the taxa table
**************************************************************************/
bool analyzeREADS::print_summary_statistics(int & num_reads)
{

    int i;
    char file[500];
    char line[1000];
    int k_value;
//    int max_index;
//    float max_val;
    float val;
    int index;
    int num_reads_dont_map = 0;
//    int mapped_reads = 0;

    ifstream in;
    ofstream out;

    int * counts=new int [num_gs];

    for(k_value=k_min; k_value<=k_max; k_value++)
    {
        for(i=0; i<num_gs; i++) counts[i]=0;
        num_reads_dont_map=0;
        sprintf(file,"%s_raw_read_results_k%d.txt",output_file_tag,k_value);
        in.open(file);
        if(!in.is_open()) {delete counts; return false;}
//        for(i=0; i<num_gs; i++) in.getline(line,300,'\t');
        in.getline(line,1000,'\n');
        num_reads=0;
        while(in.peek()!=EOF)
        {
            in.getline(line, 1000, '\t');
            if (in.peek()=='\n')
            {
                num_reads_dont_map++;
            }
            else
            {
                in.getline(line,1000,'\t');
                in.getline(line,1000,'\t');
//                max_index=-1; max_val=0.0;
                while(in.peek()!='\n')
                {
                    index = -1;
                    val = -1;
                    in>>index;
    //                index--;
                    in>>val;
    //                mapped_reads++;
    //                if(val>max_val) {max_val=val; max_index=index;}
                    counts[index]++;
                }
            }

            if (in.peek()!= EOF)
            {
                in.getline(line,300,'\n');      //get final CR
            }
//            if(max_index>=0) counts[max_index]++;
//            else num_reads_dont_map++;
            num_reads++;
        }
//        num_reads_dont_map = num_reads - mapped_reads;
        in.clear();
        in.close();
        sprintf(file,"%s_summary_statistics_k%d.txt", output_file_tag, k_value);
        out.open(file);
        out << "Number of reads:\t" << num_reads << endl;
        out << "Number of reads could not be identified:\t" << num_reads_dont_map << endl;
        out << "Taxa\tNumber of Reads" << endl;
        for(i=0; i<num_gs; i++) out << taxa[i] << "\t" << counts[i] << endl;
        out.clear();
        out.close();
    }
    delete counts;
    return true;
}

/*************************************************************************
* print_probability_statistics      //MS:broken, trouble with reading raw results
*   parameters:
*       -format_file:   0=FASTA, 1=FASTQ
*       -fof:           true=file of files (needed for formatting)
*       -file_to_open:  file name of reads
*       -k_value:       k-mer size being analyzed
*   Calculates coordinates for each read and retrieves vector of
*   probabilities based upon location in library structure l.
**************************************************************************/
bool analyzeREADS::map_file(short int format_file, bool fof, char * file_to_open, short unsigned int k_value)
{
    int i,j;
    ifstream in;
    ofstream out;
    char file[500];
    char * line=new char [5000];
    char header[1000];
    vector<lego> result;
    int factor = 100000;
    int sensi = 4;

    sprintf(file,"%s_raw_read_results_k%d.txt",output_file_tag,k_value);

    if(!fof)
    {
        out.open(file);
        out << "read" << "\t" << taxa[0];
//        out << taxa[0];
        for(i=1; i<num_gs; i++) out << "\t" << taxa[i];
//        out << endl;
    }
    else out.open(file,ios::app);

    Profile * p=new Profile(k_value,r,x_axis_file,y_axis_file);
    coordinates c;

    in.open(file_to_open);
    if (!in.is_open())
    {
        cout << "read file could not be opened" << endl;
        return false;
    }
    if(format_file==0)
    {
        while(in.peek()!=EOF)  //not finding end of file need to ask
        {
            in.getline(header,1000,'\n'); while(header[0]!='>') in.getline(header,1000,'\n');
            in.getline(line,5000,'\n');
            p->examine_read(line,c);
            result=l[k_value-k_min]->search(c, sensi);

            if(result.size()!=0)
            {
                out << endl;
                sort(result.begin(),result.end(),taxon_less());
                out << header;
                j=0;
                for(int i=0; i<num_gs; i++)
                {
                    if(result[j].taxon==i) {out << "\t" << result[j].taxon << "\t" << result[j].prob*factor; j++;}
//                    if(result[j].taxon<i) {j++; out << "0\t"; i--;}  //wat?
                }

            }
            else
            {
                out << endl << header << '\t';
                //here call to close ***TO DO***
            }
            result.clear();
            char look = in.peek();

        }
        in.clear();
        in.close();
        out.clear();
        out.close();
    }
    else
    {
        while(in.peek()!=EOF)
        {
            in.getline(header,1000,'\n'); while(header[0]!='@') in.getline(header,1000,'\n');
            in.getline(line,5000,'\n');
            p->examine_read(line,c);
            result=l[k_value-k_min]->search(c, sensi);
            if(result.size()!=0)
            {
                out << endl;
                sort(result.begin(),result.end(),taxon_less());
                out << header;
                j=0;
                for(i=0; i<num_gs; i++)
                {
                    if(result[j].taxon==i) {out << "\t" << result[j].prob*factor; j++;}  //divide by scaling factor
                    if(result[j].taxon<i) {j++; out << "0\t"; i--;}
                }

            }
            else
            {
                out << endl << header;
                //here call to close ***TO DO***
            }
            result.clear();
            in.getline(line,5000,'\n'); //sequence identifier
            in.getline(line,5000,'\n'); //quality values
        }
        in.clear();
        in.close();
        out.clear();
        out.close();
    }
    delete line;
    delete p;
    return true;
}

/*************************************************************************
* format
*   parameters:
*       -fof:           true=file of files, false=single file of reads
*   return:
*       0=FASTA, 1=FASTQ or -1 if not either
*   Determines if it is FASTA or FASTQ and if it is a file of files
**************************************************************************/
short int analyzeREADS::format(bool & fof)
{
    short int f=-1;
    char *line=new char [1000];
    ifstream in;
    in.open(file_of_files_reads);
    in.getline(line,1000,'\n');
    in.clear();
    in.close();
    if(line[0]=='>')
    {
        f=0;       //FASTA
        fof = false;
    }
    else if(line[0]=='@')
    {
        f=1;       //FASTQ
        fof = false;
    }
    else
    {
        in.open(line);
        in.getline(line, 1000, '\n');
        if (in.is_open()==true)
        {
            if (in.is_open()==true)
            {
                if (line[0]=='>')
                {
                    f=0;
                    fof = true;
                }
                if (line[0]=='@')
                {
                    f = 1;
                    fof = true;
                }
            }
        }
    }
//
//    in.open(file_of_files_reads);
//    if(in.is_open())
//    {
//        fof=true;
//        in.getline(line,1000,'\n');
//        in.clear();
//        in.close();
//    }


    delete line;
    return f;
}

/*************************************************************************
* load_library
*   Reads in library from file, sets up l.
**************************************************************************/
bool analyzeREADS::load_library()
{
    l=new legoland * [k_max-k_min+1];

    int i,j;
    short unsigned int lib_k_min, lib_k_max;
    unsigned long long multiplier;
    char line[50000];
    ifstream in;
    in.open(reference_library_file);
    if (!in.is_open())
    {
        cout << "reference library file not found" << endl;
        return false;
    }
    in >> lib_k_min;                        //gets k_min for library
    in >> lib_k_max;                        //gets k_max for library
    in.getline(line,1000,'\n');             //gets CR after k values
    in.getline(line,1000,'\n');             //gets axes files
    in>>r;
    in>>multiplier;
    in>>num_gs;
    in.getline(line,1000,'\n');             //gets CR after number of genus/species

    for(i=0; i<k_max-k_min+1; i++) l[i]=new legoland(multiplier);

    //load taxa table
    taxa=new char * [num_gs];
    for(i=0; i<num_gs-1; i++)
    {
        in.getline(line,1000,'\t');
        taxa[i]=new char [strlen(line)+1];
        strcpy(taxa[i],line);
    }
    in.getline(line,1000,'\n');
    taxa[i]=new char [strlen(line)+1];
    strcpy(taxa[i],line);

    in.getline(line,1000,'\n');             //get initial * indicating beginning of legos

    if(k_min>lib_k_min)                     //read any lower k values not wanted
    {
        for(i=lib_k_min; i<k_min; i++)
        {
            while(in.peek()!='*')
                in.getline(line,50000,'\n');
            in.getline(line,1000,'\n');     //get * separating k values
        }
    }

    //read k's
    unsigned short t;
    float p;
//    unsigned long long index_factor;
    int x,y;
    unsigned int num_values;
    lego i_val;
    std::vector<lego> val;
    for(i=0; i<k_max-k_min+1; i++)
    {
        while(in.peek()!='*')
        {
            in>>x;
            in>>y;
            in>>num_values;
            for(j=0; j<num_values; j++)
            {
                in>>t;
                in>>p;
                i_val.taxon=t;
                i_val.prob=p;
                val.push_back(i_val);
            }
            l[i]->load_legos(x, y, val);
            val.clear();
            in.getline(line,50000,'\n');                //Get CR after entry
        }
        in.getline(line,1000,'\n');     //get * separating k values
        l[i]->printLegoland();
    }
    in.clear();
    in.close();

    return true;
}

/*************************************************************************
* get_parameters
*   Get's user parameters through series of prompts. Verifies that
*   they are set properly and calls to verify files.
**************************************************************************/
bool analyzeREADS::get_parameters()
{
    if(get_reference_library())
    {
        if(get_k_range())
        {
            if(get_read_file())
            {
                get_outputfile_name();
                return true;
            }
            else return false;
        }
        else return false;
    }
    else return false;
}

/*************************************************************************
* get_reference_library
*   Gets name and location of reference library and axes files (same ones
*   used for creation of the library).
**************************************************************************/
bool analyzeREADS::get_reference_library()
{
    int gc_1, gc_2;
    char f1[500]; char f2[500];
    reference_library_file = new char [100];

    cout << "Enter the name of the Reference Library file:\n";
    cin >> reference_library_file;

    ifstream in;

    if(file_exist(reference_library_file))
    {
        in.open(reference_library_file);
        in>>k_min; in>>k_max;
        in>>gc_1; in>> gc_2;

        sprintf(f1,"GC_axis//GC%d.ref",gc_1);
        sprintf(f2,"GC_axis//GC%d.ref",gc_2);
        if(!file_exist(f1) || !file_exist(f2))
        {
            cout << "Please check that the axes files are in the folder\nlabeled GC_axes\n";
            return false;
        }

        x_axis_file=new char [strlen(f1)]; strcpy(x_axis_file,f1);
        y_axis_file=new char [strlen(f2)]; strcpy(y_axis_file,f2);
        return true;
    }
    else
        return false;
}

/*************************************************************************
* get_outputfile_name
*   Gets the name of the output file from the user
**************************************************************************/
bool analyzeREADS::get_outputfile_name()
{
    cout << "Please enter a file name for the results.\n";
    output_file_tag=new char [300];
    cin >> output_file_tag;
    return true;
}

/*************************************************************************
* get_k_range
*   Gets the k_min and k_max values from the user
**************************************************************************/
bool analyzeREADS::get_k_range()
{
    bool loop=true;
    int temp_k_min,temp_k_max;
    char selection;
    cout << "This reference library is created to consider k_min=" << k_min << " to k_max=" << k_max << endl;

    while(loop)
    {
        cout << "You can consider any k values in this range. Are these the values of k you would like to consider? (Y/N) ";
        cin>> selection;
        if(selection<96) selection+=32;

        switch(selection)
        {
            case 'y': loop=false; break;
            case 'n':
                cout<<"Enter new min k-mer size: ";
                cin>>temp_k_min;
                cout<<"Enter new max k-mer size (Maximum 9): ";
                cin>>temp_k_max;
                if( (temp_k_min<k_min) || (temp_k_max>k_max) )
                    {cout << "The library does not include these k values. Please either select new values or Q to quit.\n"; break;}
                if(temp_k_min>temp_k_max)
                    {cout << "The minimum must be less than or equal to the maximum. Please select new values.\n"; break;}
                loop=false;
                break;
            case 'q': return false;
            default: cout << "Please enter a valid selection.\n"; break;
        }
    }
    return true;
}

/*************************************************************************
* get_read_file
*   Gets the name of the read file from the user
**************************************************************************/
bool analyzeREADS::get_read_file()
{
    bool loop=true;
    cout<<"Enter the read size: ";
    cin>>r;
    file_of_files_reads = new char [100];

    while(loop)
    {
        file_of_files_reads[0] = '\0';
        cout << "Enter either a file name listing the files containing reads OR a file containing reads:\n";
        cin>>file_of_files_reads;
        if( (file_of_files_reads[0]=='q' || file_of_files_reads[0]=='Q') && (strlen(file_of_files_reads)==1) )
            return false;
        if(file_exist(file_of_files_reads)) loop=false;
        else
            cout << "That file does not exist. Please check for a valid file or enter Q to quit.";
    }
    return true;
}

/*************************************************************************
* lookup - convert the genus/species code into a genus/species name
*   parameters:
*       -code: genus/species code representation
*       -name: genus/species name as a cstring (passed by reference)
*   return: true if code found, false otherwise
**************************************************************************/
bool analyzeREADS::lookup(unsigned short code, char *& name)
{
    if(taxa[code] != NULL)
    {
        strcpy(name,taxa[code]);
        return true;
    }
    return false;
}//end lookup


////////////////////////////////////LOCAL VARIABLES////////////////////////////////////

/*************************************************************************
* file_exist - determines if a file exists
*   parameters:
*       -file_name: name of file you are checking if it exists
*   return: true if file exists else false
**************************************************************************/
bool file_exist(char * file_name)
{
    ifstream in;
    in.open(file_name);
    if(in.is_open())
    {
        in.clear();
        in.close();
        return true;
    }
    else
        return false;
}
